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COMPUTER SIMULATION OF THE EVOLUTION . OF SPIRAL GALAXIES 

By Frank Hohl 

NASA Langley Research Center 
ABSTRACT 

A two-dimensional model is used to perform computer experiments for a 
collisonless self-gravitating system. The gravitational field is obtained by 
solving the Poisson equation and the system is advanced stepwise in time. Com- 
puter simulations have been performed for systems with an initially uniform dis- 
tribution over a circular or elliptical region in x-y space, zero thermal 
velocity, and with given values of initial solid-body rotation. Up to 
10 000 stars are used in the calculations. As such systems evolve in time 
they display the various shapes observed in actual galaxies. It is found that 
some of the calculations recently reported cannot be considered collisionless. 

INTRODUCTION 

Actual stellar systems, such as galaxies, contain about 10^ stars of 
sufficiently large average separation so that binary encounters between stars 
can be neglected. Stellar systems can therefore be described by the collision- 
less Boltzmann equation. In order to simulate the evolution of such stellar 
systems on a computer one would like to follow the motion of at least several 
thousand masses or stars. Hohl and Feix (refs. 1 and 2 ) and Lecar (ref. 3 ) 
used one-dimensional sheet models to study the evolution of self-gravitating 
systems. A similar model where the motion of a large number of concentric 
spherical mass shells is followed has been used by Henon (refs. 4 and 5) • 
Recently Hockney (ref. 6) and Hohl (refs. 7 and 8) introduced a two-dimensional 
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model where the stars are represented by infinitely long mass rods. In the 
present paper the results of some simple experiments with the two-dimensional 
model are given. The calculations show that spiral and other structures similar 
to those of some actual stellar systems can be obtained by purely gravitational 
effects. 


DESCRIPTION OF THE MODEL 

The model consists of a large number of rods of mass m per unit length 
and the rods are of infinite extent in the z-direction. These rods move in 
the x-y plane under the action of their mutual gravitational attraction. 

The system of mass rods is advanced in time in the following manner. First, 
the distribution of mass p(x,y) is used to obtain the gravitational potential 
cp(x,y) by numerically solving the Poisson equation. Second, the gravitational 
field at the position of the particles is computed from the potential cp(x,y). 
Third, using Newton's laws the motion of all the mass rods is advanced for a 
constant time step 5t. These three steps represent one cycle and they are 
repeated until the desired evolution of the system is achieved. 

The crucial point in the computations is the solution of the Poisson 
equation. It is desirable that the time required for this process be only a 
small fraction of the cycle time. If the system is advanced for a small time 
step 5t the mass distribution p(x,y) will not change very much. The 
change in the gravitational potential will then also be very small. Thus, the 
solution of the finite difference form of the Poisson equation by an iteration 
method which uses the potential from the previous cycle as an initial guess 
will converge very rapidly. The accuracy of the iterative solution of the 
Poisson equation is continually checked during the calculations. This is done 
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by obtaining the values of the potential and the field at several selected 
points by summing directly the contribution from each star. The values so 
obtained agree with those obtained from the solution of the Poisson equation 
to at least the first three digits. For a 10 000 star system with a 
101 x 101 grid to solve the Poisson equation the cycle time is 7 seconds on the 
Langley Research Center's CDC 6600 data processing system. This time includes 
such operations as the checking of the gravitational potential, the calculation 
of the energy and angular momentum, and writing the positions and velocities 
of all stars on tape. 

To solve the Poisson equation one first requires the boundary conditions 
around the rectangular mesh. The potential at an arbitrary boundary point 
z = x + iy can be written as 

cp(z) = 2G £ P n ,m ln | z “ z n,m| ^ 

n,m 


where z n ^ m = x n ^ m + iy n ^ m is the coordinate of the cell n,m of the rectang- 
ular mesh, p is the mass density in that cell and G is the gravitational 
ii y m 

constant. Equation (l) can be written as 


Cp(x,y) = 2GNm ln|z| + 2G ^ ^ n>m 

n,m 


I P n,. 


Re 


In n 

z n,m\ 

V 

z / 


( 2 ) 


where N is the total number of stars in the system. Since P n ^ m is nonzero 


only for 


J n,m 


< 1, equation (2) can be written as 
cp(x,y) = 2GNm ln|z| - 2G Re 
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a k 

kz 
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( 3 ) 


3 



where 


p 

K n,m n,m 


and the series expansion for ln(l - — ^ — I has been truncated after 15 terms. 


The Poisson equation 


= taGpfx,; 


ax 2 a/ 

is solved by using the standard five-point difference equation 

^n+l,m + ^n,m+l + ^n-l^m + ^n^m-l “ ^n^m = ^ rtG Pn,m 

where Ax = Zty = 1 so that represents the number of mass rods in the 

ii y m 

cell n,m. This set of simultaneous equations is solved by an iteration method 
on the CDC 66 00 data processing system in the form 


r+1 r / r+1 r r+1 r . r \ „ \ 

~ ,m = \m + ^(Vl,m + %,l,m + Vm-1 + 9n,nH-l “ Hi,m - ^ G Pn,mj 


The superscript r refers to the r^h iteration and the parameter 7 is 
adjusted to give the maximum rate of convergence. 

The differential equations of motion for the stars are 


;(Rj,t) 


-i = V. 

dt J 


where R is the position of a star and K is the gravitational field. The 
evolution in time of the particle trajectories is given by the finite difference 


k 


app r oxima t i on 
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Rj(t + 5t) = Rj(t) + Vj(t)5t + 


V,(t + St) = V,(t) + ]?(]*,, t\&t + 

3 J \ 0 I 2 dt 


where — was approximated by the backward difference 
dt 


f ' 


K(R.,t - 5t^ 

v j y 


Equations which were found to be as accurate as equations (10) and (ll) 
but which require fewer computer operations are 


«,(* + f) - - f) + "%•* 


Hj(t + St) = Rj(t) + StV^t + 


The term — in equation (ll) was introduced to remove a computational insta- 
dt 

bility. Similarly, the use of the new velocity in equation (14) removes an 
unconditional computational instability. 

The kinetic energy of the system is 


T - 1 m L h 2 


and the potential energy is 


IN 

- 1 " Z <{®j) 


where the summation goes to N, the number of stars in the system. The total 
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energy of the system U remains constant and is 


U = T + p (17) 

The relaxation time (ref. 9) for the two-dimensional model is of interest 
to determine for how long the system can he treated as collisionless. 

Consider an encounter between a rod of mass m and velocity V with a 
stationary rod of mass M. If D is the distance of closest approach (impact 
parameter) then the transverse force acting on the moving star during a time 
2D/V can be approximated as 2GfflM/D. The moving star has therefore acquired 
a transverse velocity 


5V = 4GM/V 


(18) 


For a star interacting with a system containing many stars the effect of many 
individual encounters must he summed. The number of encounters in a time t 
with impact parameter between D and D + dD is tVn dD where n is the 

p 

density of stars in the system. The expectation value of is then given 

fey 


<Vj_ 2 > = 



l6tnG%^ 

V 


dD 


- l6tnO%l2(D nax - D mln )/v 

- l6tnG 2 M 2 D max /V ( 19 ) 


since in general D max » ^min where D max is the dimension of the system. 

o 

The relaxation time r Q is defined to be the time required for y to be 

of the same order of magnitude as V^. The velocity V is taken to be the 
characteristic velocity V g = h max / T g where T g = 2* // 2rtGmn is the rotation 
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period of a uniform circular distribution of mass rods with the centrifugal 
force equal to the gravitational attraction. Equation (19) then becomes 

V 5 

t c = g (20) 

l6nG 2 m\ iax 

where M = m. The ratio of the relation time to the rotation period is then 

1 

I given by 


l£ s JL 

T g 500 


( 21 ) 


Equation (21) is only a rough estimate but it nevertheless indicates that for 
N equal to a few thousand, colli sional effects can be expected to become 
important after only a few rotations. 


RESULTS AND DISCUSSION 

Computer simulations have been performed for a number of systems with an 
initially uniform distribution over a circular or elliptical region in x-y 
coordinate space, zero thermal velocity, and various values of initial solid- 
body rotation. In particular, systems with a ratio of minor axis a to major 
axis b equal to 3, and with an initial solid-body rotation given by 

(JOg 2 = 8NmG /(a + b) 2 (22) 

have been investigated in detail. It was found that for a system containing 
2000 stars with a 51 X 51 grid to solve the Poisson equation, eollisional 
effects caused sufficient randomization to destroy the spiral structure after 
only 2 to 3 rotations. The results are shown in figure 1 to 3* The normaliza- 
tion G = m = 1 were used in the calculations and the time is shown in rota- 
tion periods t_ = 2rt /<£>„ . Figure 1 shows the evolution of the 2000-star system 
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in x-y coordinate space. The evolution of the velocity distribution in 
V x - Vy space is shown in figure 2. Figure 5 shows the evolution of the 
velocity distribution after the initial solid body rotation has subtracted. 

The radial velocity V r is plotted versus ^Vq - ro)g^ where V Q is the azi- 
muthal velocity and r is the radius from the center of the system to a star. 
Initially the thermal (or random) velocity of all mass rods is zero. As the 
system evolves in time the thermal velocity builds up and the system expands in 
v r - (V e - roJj space. After a time t = 1.75Tg there is little further change 
in the velocity distribution. Thus, the calculations reported by Hockney 
(ref. 6) with 2000 stars and a 48 x 48 grid cannot be considered collisionless 
for more than 2 rotations. The same holds for the results given by Hohl (ref. 8) 
for 2000 stars (with a 51 X 51 grid). To reduce collisional effects the calcu- 
lations were repeated with an identical initial elliptical distribution but 
with 4000 stars and a 101 x 101 grid to solve the Poisson equation. The results 
are shown in figures 4 and 5- Figure 4 shows that the spiral structure now has 
a somewhat greater stability but it still tends to randomize after about 
4 rotations. The corresponding evolution in V x - Vy velocity space is shown 
in figure 5* The number of stars was increased to 10,000 and the resulting 
evolution in x-y coordinate space is shown in figure 6. The evolution of yet 
another system with 20,000 stars was computed up to four rotations and it 
resulted in a structure nearly identical to that shown in figure 6. Thus, it 
appears that in order to simulate collisionless systems for the first n rota- 
tions, nearly 2n X 10^ stars are required. Figure 1 shows that at t = 6Tg 
the spiral structure begins to randomize. To keep the system collisionless for 
a longer time calculations are now being performed with systems containing up 
to 200,000 stars. 
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An interesting point is that all systems investigated which have an 
initially elongated distribution and a given solid-body rotation go through an 
initial evolution similar to that shown in figure 6 up to t = I.Ot . At 

o 

t = 0.75T g the system shows a structure nearly identical to that of the 
peculiar galaxy (GBl) discussed by Burbidge, Burbidge, and Shelton (ref. 10). 
Another such integral-sign shaped galaxy is NGC 4656/7 reproduced in the Hubble 
Atlas (ref. 11). Figure 6 shows that at t = 3 . 0 t the two spiral arms begin 

O 

to overlap and give the appearance of a circular ring separate from the 
nucleus, presenting a shape similar to NGC 514-5 (ref. 11). At t = 4.0r g the 
system shows a structure similar to NGC 5194 (ref. 11). An interesting point 
is that at t = 3*C>T g the spiral structure has almost disappeared and at 
t = 3*75 T gj after less than one rotation, the spiral structure is again very 
pronounced. This phenomena appears to reoccur from t = 4.5C>T g to 
t = 6.00r g . However, at t = 6.00T g collisional effects (or computer noise) 
have already appreciably affected the evolution of the system. The author 
expects that for the calculations presently in progress with 100,000 stars, 
■collisional effects will be small enough so that several such transitions from 
a ring shape to a spiral structure can be observed. Such repeating density 
waves resulting in spiral structure will prevent the arms from winding up. 

This overall picture is in agreement with the theory of Lin and Shu (ref. 12) 
where the stars pile up for a certain time in the spiral-shaped potential 
wells. The evolution of the velocity distribution in V x - Vy space is 
shown in figure J. At t = 0.25T g the system shows an interesting structure. 
Condensation in velocity space has occurred and striation appears in the 
velocity distribution. The corresponding evolution in V r - (Vq - roA 
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velocity space is shown in figure 8. This distribution shows little change 
after the first 3 rotations. 

One might expect that the appearance of the two-arm spiral structure 
shown in figure 6 depends on the initially elongated distribution. This is 
not the case as is illustrated in figure 9 « Figure 9 shows the evolution of 
a 10,000-star system with an initially uniform distribution over a circular 
region and with an initial solid body rotation equal to that required to 
balance gravitation, that is, cu = ojg = UrtGNm. For the first two rotations 
the evolution of the balanced cylinder is very similar to that of the 
2000-star system presented by Hockney (ref. 6) and Hohl (ref. 8). Figure 10 
shows the evolution of such a 2000-star system in x-y coordinate space. 

For the 2000- star system it is found that the relatively small fourth mode 
surface wave around the periphery of the cylinder quickly change to a second 
mode surface wave (egg shape), then all structure is lost and the system 
finally assumes a circular shape. Figure 9 shows that the evolution is quite 
different when colli sional effects are reduced. The fourth mode surface wave 
now continues to grow and after five rotations a four-arm spiral begins to 
form. After about 8 rotations two of the arms have disappeared and the system 
has assumed a bar shape. One would now expect that a two-arm spiral structure 
would form. However, figures 11 and 12 shows that after about nine rotations 
colli sional effects have increased the random or thermal velocity to such an 
extent that they are of the same magnitude as the rotational velocities and the 
model no longer simulates a collisionless system. Nevertheless, figure 9 shows 
that the second mode surface wave will finally dominate. Figure 11 shows that 
there is little change of the distribution in V x - Vy velocity space for the 
first 5 rotations. After that time condensation in velocity space begin to 
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occur. The evolution of the random velocity shown in figure 12 shows that 
there is an almost linear increase in the random or thermal motion of the mass 
rods. 

It would be more appropriate to simulate spiral galaxies by means of a 
model using finite length mass rods or even point masses moving in a plane. 

The two-dimensional model does simulate spiral structure and shows that two-arm 
spirals are very likely to occur. However, whether the two-dimensional model 
correctly simulates some of the more refined properties of spiral galaxies can 
only be determined by repeating the calculations with a model of finite length 
rods or point masses. The reason for using the two-dimensional model is that 
the solution of the Poisson equation is relatively simple. Since it was found 
that the number of stars required to simulate a collisonless system for a suf- 
ficiently long time is of the order of 100 000, the time used for solving the 
Poisson equation represents only a very small portion of the cycle time. The 
additional time required for obtaining the potential distribution for point 
masses will, therefore, no longer be an important factor. For this reason 
future calculations will be performed for systems of finite length rods and 
point masses moving in a plane. 
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Figure 1. - Evolution of a 2000-star system in coordinate space. 
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(a) Evolution up to t = 2.75^g. 


Figure 4, - Evolution of a 4000-star system in x-y coordinate space 
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(a) Evolution up to t = 2.75"Tg. 

Figure 5- - Evolution of a 4000-star system in V-*- - V v velocity space 
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(a) Evolution up to t = 2.75Tg* 

Figure 9«- Evolution of a 10, 000-star system in x-y coordinate space. 
The initial solid body rotation of the system equalled that required 
to balance the gravitational attraction. 


















Figure 10.- Evolution of a 2,000-star system where initially the 
centrifugal force and gravitational attraction are in balance 



t = 0 t = 0.25 T t = 0.50 r 

s g 



t = 2.25 t t = 2.50 T t = 2.75 r 

S g g 


(a) Evolution up to t = 2.75Tg. 

Figure 11.- Evolution of a 10j000-star system in velocity wher 
initially the centrifugal force balances the gravitational 
attraction. 
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(a) Evolution up to t = 2.75Tg. 

Figure 12.- Evolution of a 10,000-star system in V r - ^Vq - vel 

space. The initial solid body rotation of the system equalled that 
required to balance the gravitational attraction. 










